#ifndef OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
#define OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_center.h>
#include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_minimum_coord.h>
#include <OpenTissue/core/containers/mesh/common/util/mesh_compute_mesh_maximum_coord.h>
#include <OpenTissue/core/containers/mesh/common/util/mesh_deformation_modifiers.h>

#include <cassert>
#include <vector>

namespace OpenTissue
{
  namespace mesh
  {

    /*
    First, I will assume you are more or less familiar with the Marching Cubes algorithm.

    The cube description I use looks like this:

    Y
    |
    3----2----2
    /|        /|
    11 |      10 |
    /  3      /  1
    7----6----6   |
    |   |     |   |
    |   0----0|---1---X
    7  /      5  /
    | 8       | 9
    |/        |/
    4----4----5
    /
    Z

    Given that, here's the basic procedure I have developed:

    You know that Marching Cubes treats a volume as a regular 3D grid of "cubes", with
    8 voxel values of the volume at the corners. For each of the 12 edges in any given cube,
    the isovalue can fall above, below, or between the values at the endpoints. This means
    that any edge can have at most one interpolated vertex "on" it.

    The main data structure for my extraction algorithm is a "Cell" illustrated by the
    following diagram:

    Y
    |
    |
    Vy
    |
    |
    Vs---Vx----X
    /
    Vz
    /
    Z

    The "CellInfo" struct below represents this structure. If you imagine putting lots
    of these cells next to each other in space, you can quickly see they form a grid similar
    to a 3D grid of cubes, but with no duplicated edges:

    |   |   |
    +-- +-- +--
    /   /   /
    |   |   |
    +-- +-- +--
    /   /   /

    Two 2D arrays of cells, each the same size as a "slice" in the 3D volume, will between them contain
    all necessary information for any vertices and triangles generated by the first set of cubes
    in the volume. Each 2D array of cells I'll from hereon refer to as a "grid".

    Looking at the Cell diagram above and CellInfo struct below:
    Vx is the integer index into a list of vertices for the vertex interpolated along that edge - the x-axis.
    That is, if the isovalue falls between the values of the endpoints of that edge, then a vertex is
    interpolated and stored in a list, and its index in that list is stored in Vx.
    It is -1 initially, or if no vertex has yet been interpolated. Likewise for the other 2 axes.

    During extraction, each "cube" of 8 values are visited, and used to construct the index into
    'triangleEdgeTable', which has the list of triplets of edges forming at most 5 triangles for
    the cube. As each of the 3 edges forming the triangle are visited, the "EdgeCellInfo" struct
    is used to see which cell that edge maps to in one of the 2 grids.

    For example, edge 9 in one cell maps to Vz in the cell next to it, in the same grid; i.e. to
    the right of it in the cube diagram above. Edge 5 in one cell maps to Vy in the cell to the
    right of it and in the next grid.

    Now that we know which cell the edge maps to, we can see if the vertex on that edge has been
    interpolated by checking to see if the corresponding 'vertIndex' field is -1. If it is, we
    interpolate the vertex, add it to the list, and store its index. If not, then the vertex has
    already been interpolated and we just use it.

    We now have the 3 indexes of the vertices that form a triangle, and can add that triangle to the
    list of triangles. The normal of the triangle is computed and stored as well as accumulated
    within each vertex. Later, the per-vertex normals are averaged and normalized.

    After processing 2 slices of the volume, the "top" grid (i.e. the 2nd grid) can be reused as the
    "bottom" grid for the next set of slices. It is swapped with the 1st grid, and the new 2nd grid is
    reset and reused.

    I have also implemented a vertex "smoothing" operation, which basically collects vertices that
    are within 1/2 unit of a "corner" of the cube and averages them, forming as few as half as many
    vertices and triangles and preserving the overall shape and morphology of the surface.

    The 'smoothedIndex' field holds a pointer to the smoothed vertex in the list, and the 'vertIndex'
    holds pointers to which smoothed vertex is closest to the corresponding edge. Everything else
    works the same, except degenerate triangles (triangles where two or more of the vertices are
    the same) are not added to the triangle list.

    At the end, smoothed vertices are averaged, and all vertices are pre-scaled according to a caller-specified
    scaling factor (to account for isotropism) and translated so as to be in a default [-0.5,0.5] space
    centered at the origin.

    */

    namespace detail
    {

      /* Marching Cubes Isosurface extraction routines.
      The Marching Cubes algorithm is patented and was originally developed
      by William Lorenson and Harvey Cline:

      "Marching Cubes: A High Resolution 3D Surface Construction Algorithm",
      William E. Lorensen and Harvey E. Cline,
      Computer Graphics (Proceedings of SIGGRAPH '87), Vol. 21, No. 4, pp. 163-169.
      */

      /* The triangleEdgeTable was obtained from open-source code written by Cory Bloyd
      which was based on a marching cubes description by Paul Bourke:
      http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/

      Everything else in here is mine! :-)
      */



      /**
      * this holds per-voxel info during isosurface extraction.
      *  2 2D grids are allocated the same size as a 2D "slice" of the volume
      */
      template<typename mesh_type>
      class CellInfo
      {
      public:

        typedef typename mesh_type::vertex_handle vertex_handle;

      public:

        vertex_handle  m_vertex_handle[3];  ///< The handle  of the vertex interpolated on each axis (Vx,Vy,Vz)
        vertex_handle  m_smoothed_handle;   ///< The handle of the smoothed (accumulated) vertex

      public:

        void clear()
        {
          vertex_handle null_handle;
          m_vertex_handle[0]  = null_handle;
          m_vertex_handle[1]  = null_handle;
          m_vertex_handle[2]  = null_handle;
          m_smoothed_handle = null_handle;
        }
      };

      /**
      * for each edge in a cube, this tells which axis in which cell in which grid the edge maps to,
      *  and its spacial coordinate offset
      */
      class EdgeCellInfo
      {
      public:
        int m_cell;
        int m_grid;
        int m_axis;
        int m_offset[3];
      };

      template<typename cell_info_iterator>
      void clear_cell_grid(cell_info_iterator begin,cell_info_iterator end)
      {
        for(cell_info_iterator c=begin;c!=end;++c)
          c->clear();
      }

    }// namespace detail

    /**
    * Extract Iso Surface.
    *
    * @param phi
    * @param isolevel
    * @param mesh
    * @param sub_sample
    * @param smoothing
    */
    template<typename grid_type,typename mesh_type>
    void smooth_isosurface(
      grid_type const & phi
      , typename grid_type::value_type const & isolevel
      , mesh_type & mesh
      , size_t sub_sample  = 1u
      , bool smoothing   = true
      )
    {
      typedef detail::CellInfo<mesh_type>           cell_info_type;
      typedef std::vector<cell_info_type>           cell_grid_type;
      typedef typename grid_type::value_type         value_type;

      typedef typename mesh_type::math_types        math_types;
      typedef typename math_types::value_traits     value_traits;
      typedef typename math_types::vector3_type     vector3_type;
      typedef typename math_types::real_type        real_type;

      typedef typename mesh_type::vertex_handle     vertex_handle;
      typedef typename mesh_type::vertex_iterator   vertex_iterator;

      static int triangle_edge_table[256][16] =
      {
        {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
        {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
        {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
        {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
        {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
        {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
        {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
        {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
        {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
        {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
        {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
        {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
        {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
        {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
        {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
        {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
        {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
        {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
        {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
        {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
        {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
        {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
        {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
        {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
        {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
        {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
        {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
        {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
        {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
        {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
        {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
        {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
        {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
        {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
        {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
        {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
        {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
        {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
        {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
        {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
        {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
        {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
        {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
        {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
        {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
        {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
        {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
        {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
        {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
        {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
        {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
        {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
        {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
        {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
        {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
        {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
        {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
        {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
        {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
        {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
        {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
        {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
        {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
        {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
        {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
        {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
        {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
        {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
        {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
        {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
        {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
        {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
        {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
        {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
        {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
        {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
        {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
        {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
        {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
        {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
        {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
        {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
        {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
        {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
        {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
        {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
        {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
        {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
        {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
        {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
        {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
        {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
        {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
        {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
        {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
        {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
        {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
        {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
        {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
        {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
        {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
        {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
        {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
        {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
        {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
        {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
        {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
        {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
        {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
        {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
        {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
        {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
        {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
        {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
        {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
        {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
        {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
        {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
        {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
        {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
        {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
        {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
        {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
        {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
        {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
        {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
        {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
        {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
        {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
        {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
        {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
        {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
        {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
        {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
        {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
        {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
        {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
        {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
        {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
        {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
        {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
        {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
        {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
        {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
        {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
        {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
        {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
        {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
        {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
        {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
        {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
        {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
        {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
        {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
        {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
        {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
        {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
        {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
        {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
        {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
        {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
        {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
        {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
        {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
        {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
        {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
        {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
        {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
        {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
        {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
        {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
        {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
        {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
        {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
        {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
        {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
        {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
        {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
        {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
        {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
        {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
      };


      /**
      * which vertices in the cube make up each edge?
      */
      static int edge_vertex_table[12][2] =
      {
        {0,1}, {1,2}, {3,2}, {0,3},
        {4,5}, {5,6}, {7,6}, {4,7},
        {0,4}, {1,5}, {2,6}, {3,7}
      };

      //--- Clear any data stored in mesh.
      mesh.clear();

      //--- Determine constants used by the algorithm
      if( sub_sample <= 0 )  sub_sample = 1;
      if( sub_sample > 16 )  sub_sample = 16;
      size_t  sampled_I   = phi.I()/sub_sample;
      size_t  sampled_J   = phi.J()/sub_sample;
      //size_t  sampled_K   = phi.K()/sub_sample;
      size_t  grid_width  = sampled_I-1;
      size_t  grid_height = sampled_J-1;
      // NOTE: What about some sort of warning here for the user? An assertion perhaps.
      if( grid_width == 0 || grid_height == 0 )
        return;

      size_t cnt_cells = sampled_I * sampled_J;

      //--- Determine offsets needed by the algorithm...
      detail::EdgeCellInfo     edge_cell_lut[12];
      {
        //--- compute the edge-cell lookup: given any edge in the cube (there are 12),
        //--- this table indicates which cell the interpolated vertex for that edge is
        //--- (as indexed from the current cell), which grid that cell is in (either the
        //--- first or the second), and which axis the edge maps to (either x, y, or z).
        //--- grid is which grid of CellInfos to use: the first or the second
        //--  cell is how to index from the current cell to the proper cell for that edge
        //--- axis is which axis the needed vertex is on: 0=>x, 1=>y, 2=>z
        edge_cell_lut[ 0].m_cell = 0;          edge_cell_lut[ 0].m_grid = 0;    edge_cell_lut[ 0].m_axis = 0;
        edge_cell_lut[ 1].m_cell = 1;          edge_cell_lut[ 1].m_grid = 0;    edge_cell_lut[ 1].m_axis = 1;
        edge_cell_lut[ 2].m_cell = sampled_I;   edge_cell_lut[ 2].m_grid = 0;    edge_cell_lut[ 2].m_axis = 0;
        edge_cell_lut[ 3].m_cell = 0;          edge_cell_lut[ 3].m_grid = 0;    edge_cell_lut[ 3].m_axis = 1;
        edge_cell_lut[ 4].m_cell = 0;          edge_cell_lut[ 4].m_grid = 1;    edge_cell_lut[ 4].m_axis = 0;
        edge_cell_lut[ 5].m_cell = 1;          edge_cell_lut[ 5].m_grid = 1;    edge_cell_lut[ 5].m_axis = 1;
        edge_cell_lut[ 6].m_cell = sampled_I;   edge_cell_lut[ 6].m_grid = 1;    edge_cell_lut[ 6].m_axis = 0;
        edge_cell_lut[ 7].m_cell = 0;          edge_cell_lut[ 7].m_grid = 1;    edge_cell_lut[ 7].m_axis = 1;
        edge_cell_lut[ 8].m_cell = 0;          edge_cell_lut[ 8].m_grid = 0;    edge_cell_lut[ 8].m_axis = 2;
        edge_cell_lut[ 9].m_cell = 1;          edge_cell_lut[ 9].m_grid = 0;    edge_cell_lut[ 9].m_axis = 2;
        edge_cell_lut[10].m_cell = 1+sampled_I; edge_cell_lut[10].m_grid = 0;    edge_cell_lut[10].m_axis = 2;
        edge_cell_lut[11].m_cell = sampled_I;   edge_cell_lut[11].m_grid = 0;    edge_cell_lut[11].m_axis = 2;

        //--- The offset is how much to add to the coordinates of the current
        //--- cell to get the right interpolated  vertex coordinate
        edge_cell_lut[ 0].m_offset[0] = 0; edge_cell_lut[ 0].m_offset[1] = 0; edge_cell_lut[ 0].m_offset[2] = 0;
        edge_cell_lut[ 1].m_offset[0] = 1; edge_cell_lut[ 1].m_offset[1] = 0; edge_cell_lut[ 1].m_offset[2] = 0;
        edge_cell_lut[ 2].m_offset[0] = 0; edge_cell_lut[ 2].m_offset[1] = 1; edge_cell_lut[ 2].m_offset[2] = 0;
        edge_cell_lut[ 3].m_offset[0] = 0; edge_cell_lut[ 3].m_offset[1] = 0; edge_cell_lut[ 3].m_offset[2] = 0;
        edge_cell_lut[ 4].m_offset[0] = 0; edge_cell_lut[ 4].m_offset[1] = 0; edge_cell_lut[ 4].m_offset[2] = 1;
        edge_cell_lut[ 5].m_offset[0] = 1; edge_cell_lut[ 5].m_offset[1] = 0; edge_cell_lut[ 5].m_offset[2] = 1;
        edge_cell_lut[ 6].m_offset[0] = 0; edge_cell_lut[ 6].m_offset[1] = 1; edge_cell_lut[ 6].m_offset[2] = 1;
        edge_cell_lut[ 7].m_offset[0] = 0; edge_cell_lut[ 7].m_offset[1] = 0; edge_cell_lut[ 7].m_offset[2] = 1;
        edge_cell_lut[ 8].m_offset[0] = 0; edge_cell_lut[ 8].m_offset[1] = 0; edge_cell_lut[ 8].m_offset[2] = 0;
        edge_cell_lut[ 9].m_offset[0] = 1; edge_cell_lut[ 9].m_offset[1] = 0; edge_cell_lut[ 9].m_offset[2] = 0;
        edge_cell_lut[10].m_offset[0] = 1; edge_cell_lut[10].m_offset[1] = 1; edge_cell_lut[10].m_offset[2] = 0;
        edge_cell_lut[11].m_offset[0] = 0; edge_cell_lut[11].m_offset[1] = 1; edge_cell_lut[11].m_offset[2] = 0;
      }
      //--- Notice the convention that lowercase i,j,k are local subsampled
      //--- indices, uppercasre I,J,K are the non-sampled indices
      //--- (ie. the ones you should use to access the value phi(i,j,k))
      //--- The conversion formulas are:
      //---
      //---  I = i*sub_sample
      //---  J = j*sub_sample
      //---  K = k*sub_sample
      //---
      //--- Compute offsets to access values of nodes in phi. That is:
      //---
      //---  phi(           I,           J,K)
      //---  phi(I+sub_sample,           J,K)
      //---  phi(I+sub_sample,J+sub_sample,K)
      //---  phi(           I,J+sub_sample,K)
      //---
      size_t offset[4];
      {
        offset[0] = 0;                       //--- offset to subsampled location (  i,j,k)
        offset[1] = sub_sample;              //--- offset to subsampled location (i+1,j,k)
        offset[2] = (1+phi.I())*sub_sample;  //--- offset to subsampled location (i+1,j+1,k)
        offset[3] = phi.I()*sub_sample;      //--- offset to subsampled location (  i,j+1,k)
      }

      //--- Allocate memory for internal data structures and initialize any internal data
      static cell_grid_type   grid_storage[2];     //--- The actual memory used to store the grid cells!!!
      cell_info_type * cell_grid[2];        //--- Used to swap cell grids when moving on to next z-plane.
      cell_info_type * grid[2];             //--- Used to traverse cells in the grids!!!
      value_type  const * data = phi.data();   //--- pointer to location phi(0,0,0)
      {
        grid_storage[0].resize( cnt_cells );
        grid_storage[1].resize( cnt_cells );
        cell_grid[0] = &(grid_storage[0][0]);
        cell_grid[1] = &(grid_storage[1][0]);
        clear_cell_grid(cell_grid[0],cell_grid[0]+cnt_cells);
      }

      size_t k,K;
      for( k = 0, K = 0; K < (phi.K()-sub_sample) ; ++k, K+=sub_sample ) //--- run over sub-sampled z-planes
      {
        grid[0] = cell_grid[0];
        grid[1] = cell_grid[1];
        clear_cell_grid(cell_grid[1],cell_grid[1]+cnt_cells);
        for(size_t j = 0; j < grid_height; ++j )
        {
          value_type const * vox0 = data + j*sub_sample*phi.I() + K*phi.I()*phi.J();                //--- pointer to sub sampled location phi(i,j,k),    with i=0
          value_type const * vox1 = data + j*sub_sample*phi.I() + (K+sub_sample)*phi.I()*phi.J();   //--- pointer to sub sampled location phi(i,j,k+1),  with i=0

          for(size_t i = 0; i < grid_width; ++i )
          {
            //--- Compute the index for which set of triangles to use
            int index = 0;
            if( (float)(vox0[offset[0]]) <= isolevel )      index |= 1;    //--- phi(  i,  j,k  )
            if( (float)(vox0[offset[1]]) <= isolevel )      index |= 2;    //--- phi(i+1,  j,k  )
            if( (float)(vox0[offset[2]]) <= isolevel )      index |= 4;    //--- phi(i+1,j+1,k  )
            if( (float)(vox0[offset[3]]) <= isolevel )      index |= 8;    //--- phi(  i,j+1,k  )
            if( (float)(vox1[offset[0]]) <= isolevel )      index |= 16;   //--- phi(  i,  j,k+1)
            if( (float)(vox1[offset[1]]) <= isolevel )      index |= 32;   //--- phi(i+1,  j,k+1)
            if( (float)(vox1[offset[2]]) <= isolevel )      index |= 64;   //--- phi(i+1,j+1,k+1)
            if( (float)(vox1[offset[3]]) <= isolevel )      index |= 128;  //--- phi(  i,j+1,k+1)
            //--- Quick check for the most common indexes of all (i.e. no triangles formed)
            if( index==0 || index==255 )
            {
              //--- Advance to next pixel/cell
              vox0 += sub_sample;
              vox1 += sub_sample;
              grid[0]++;
              grid[1]++;
              continue;
            }
            int * triangle_edge = triangle_edge_table[index];
            for(int t = 0; ; t += 3 )
            {
              if( triangle_edge[t] == - 1 )
                break;
              std::vector<vertex_handle> vertex_handle(3);
              //--- Examine each of the three edges of the cell
              for( int e = 0; e < 3; ++e )
              {
                int E = triangle_edge[t+e];                                   //--- E???
                detail::EdgeCellInfo * edge = edge_cell_lut + E;              //--- this is just a faster version of: edge = &(edge_cell_lut[E]);
                int axis = edge->m_axis;                                      //--- the axis on where the intersection lies
                cell_info_type * cell = &(grid[edge->m_grid][edge->m_cell]);  //--- get cell corresponding to the intersected edge
                if( cell->m_vertex_handle[axis].is_null() )                   //--- see if we already got a vertex on this edge
                {
                  //--- Get the 2 vertices that make up the edge.
                  //--- These values are on [0,7], where values [0,3] refer
                  //--- to grid 0, and [4,7] refer to grid 1
                  int v0 = edge_vertex_table[E][0];
                  int v1 = edge_vertex_table[E][1];
                  //--- get the corresponding phi values of the 2 vertices
                  value_type val0 = ( v0 < 4 ) ? vox0[offset[v0]] : vox1[offset[v0-4]];
                  value_type val1 = ( v1 < 4 ) ? vox0[offset[v1]] : vox1[offset[v1-4]];
                  //--- Compute the interpolant (i.e. the fraction the new vertex is between the 2 others)
                  value_type frac = (isolevel-val0)/(val1-val0);
                  //--- Add a new vertex to the list of vertices

                  vector3_type coord(0,0,0);
                  cell_info_type * new_cell = 0;

                  switch( axis )
                  {
                  case 0:
                    coord(0) = static_cast<real_type>((i + edge->m_offset[0]) + frac);
                    coord(1) = static_cast<real_type> (j + edge->m_offset[1]);
                    coord(2) = static_cast<real_type> (k + edge->m_offset[2]);
                    new_cell = &(grid[edge->m_grid][edge->m_cell+1]);
                    break;
                  case 1:
                    coord(0) = static_cast<real_type> (i + edge->m_offset[0]);
                    coord(1) = static_cast<real_type>((j + edge->m_offset[1]) + frac);
                    coord(2) = static_cast<real_type> (k + edge->m_offset[2]);
                    new_cell = &(grid[edge->m_grid][edge->m_cell+sampled_I]);
                    break;
                  default:
                    coord(0) = static_cast<real_type> (i + edge->m_offset[0]);
                    coord(1) = static_cast<real_type> (j + edge->m_offset[1]);
                    coord(2) = static_cast<real_type>((k + edge->m_offset[2]) + frac);
                    new_cell = &(grid[edge->m_grid+1][edge->m_cell]);
                    break;
                  }

                  //--- if we're doing vertex smoothing, we need to accumulate
                  //--- all potential verts into a single "smoothed" vert
                  if( smoothing )
                  {
                    if( frac <= static_cast<value_type>(0.5) )//--- hmmm, why this?
                    {
                      new_cell = cell;
                    }
                    if( new_cell->m_smoothed_handle.is_null() )
                    {
                      //--- the smoothed vertex has not been started yet. Make the interpolated vertex the first
                      new_cell->m_smoothed_handle = mesh.add_vertex( coord);
                      vertex_iterator v = mesh.get_vertex_iterator(new_cell->m_smoothed_handle);
                      v->m_tag = 1;  //--- use tag to remember the  number of vertices that have been smoothed
                      cell->m_vertex_handle[axis] = new_cell->m_smoothed_handle;
                      vertex_handle[e] = new_cell->m_smoothed_handle;
                    }
                    else
                    {
                      //--- we have already started smoothing. Accumulate the interpolated vertex
                      vertex_iterator v = mesh.get_vertex_iterator(new_cell->m_smoothed_handle);
                      v->m_tag++;               //--- increment the number of smoothed vertices...
                      v->m_coord += coord;
                      vertex_handle[e] = new_cell->m_smoothed_handle;
                      cell->m_vertex_handle[axis] = new_cell->m_smoothed_handle;
                    }
                  }
                  else
                  {
                    //--- store the index of the new vertex into the cell, and remember it for this edge
                    vertex_handle[e] = mesh.add_vertex(coord);
                    vertex_iterator v = mesh.get_vertex_iterator(vertex_handle[e]);
                    v->m_tag = 1;
                    cell->m_vertex_handle[axis] = vertex_handle[e];
                  }
                }
                else
                {
                  //--- This edge has already been interpolated, so just use it
                  vertex_handle[e] = cell->m_vertex_handle[axis];
                }
              }
              //--- Add a new triangle
              if( vertex_handle[0] == vertex_handle[1] || vertex_handle[1] == vertex_handle[2] || vertex_handle[0] == vertex_handle[2] )
              {
                //--- It is weird to not have anything in here, but degenerate triangles are discovered faster  due to lazy evaluation.
              }
              else
              {
                mesh.add_face( vertex_handle[0], vertex_handle[1], vertex_handle[2] );
              }
            }
            //--- advance to next pixel/cell
            vox0 += sub_sample;
            vox1 += sub_sample;
            grid[0]++;
            grid[1]++;
          }
          //--- Advance to next row (have to increment once more because loop stops at sampled_I-1)
          grid[0]++;
          grid[1]++;
        }
        //--- Swap the two grids as we advance to the next slice
        cell_info_type * tmp = cell_grid[0];
        cell_grid[0] = cell_grid[1];
        cell_grid[1] = tmp;
      }

      {
        vertex_iterator  v = mesh.vertex_begin();
        vertex_iterator  vend = mesh.vertex_end();
        for(;v!=vend;++v)
          if(v->m_tag>1)
          {
            real_type N = static_cast<real_type>(1.0/v->m_tag);
            v->m_coord *= N;
          }
      }

      mesh::scale(mesh,sub_sample*phi.dx(),sub_sample*phi.dy(),sub_sample*phi.dz());
      mesh::translate(mesh,phi.min_coord());

      mesh::compute_angle_weighted_vertex_normals(mesh);
    }

  } // namespace mesh
} // namespace OpenTissue

//OPENTISSUE_CORE_CONTAINERS_MESH_COMMON_UTIL_MESH_SMOOTH_ISOSURFACE_H
#endif
